%% Load Mackey-Glass timeseries data
load mgdata.dat
time = mgdata(:, 1); x = mgdata(:, 2);
%x = sin(time/50);
%figure(1), plot(time, x);
%title('Mackey-Glass Chaotic Time Series')
%xlabel('Time (sec)')
%% use krls (Wingate)
kernel_func = @rbf4nn;   % the Gaussian kernel
kparam = 100;              % the variance of the Gaussian
ald_thresh = 0.1;        % the linear dependence threshold
num_samples = 100;
num_dims = 1;
% Run the KRLS algorithm
kp = krls_init( kernel_func, kparam, ald_thresh, time(1), x(1) );
% feed in examples one at a time
tic
for i = 2:num_samples
  %disp(num2str(i));
  kp = krls(kp, time(i), x(i));
end;
toc
% predict
maxP = 100;
xStar = zeros(num_dims,maxP);
for i = 1 : maxP % size(time,1) - num_samples
    xStar(i) = krls_query(kp, time(num_samples+i));
end
figure(1), plot(time, x);
title('Mackey-Glass Chaotic Time Series')
xlabel('Time (sec)')
hold on
plot(time(num_samples+1:num_samples+maxP),xStar,'r+');
hold off
%% use krls (Engel)
% write training file
clear dic nDic
nTrain = 100;
dlmwrite('train.dat', [nTrain 2], 'delimiter', ' ')
dlmwrite('train.dat', [time(1:nTrain) x(1:nTrain)], '-append', 'roffset', 0, 'delimiter', ' ')
% do the krls
!./krls train.dat
type train.dat.out
% read the result
fid = fopen('train.dat.out');
nDic = fscanf(fid, '%g', [1])
dic = fscanf(fid, '%g %g', [nDic 2])
fclose(fid);
%% do the forecasting
% 1. compute the kernel
k = rbf4nn( (101:200)', dic(:,1), 1 )
% 2. compute the weighted sum (which is the forecast) 
kp.dp.Dict = dic(:,1)';
kp.Alpha = dic(:,2:end);
for i = 1 : maxP % size(time,1) - num_samples
    xStar(i) = krls_query(kp,  nTrain+1);
end
figure(1), plot(time, x);
title('Mackey-Glass Chaotic Time Series')
xlabel('Time (sec)')
hold on
plot(time(num_samples+1:num_samples+maxP),xStar,'r+');
hold off


